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1.  Introduction 

Phase  field  models  have  been  successfully  applied  towards 
many  problems  in  continuum  physics.  At  each  point  in  the  problem 
domain,  one  or  more  phase  field  variables  describe  the  state 
of  the  substance.  Conserved  phase  field  variables  are  typically 
related  to  composition,  e.g.,  molar  or  mass  fractions  of  atoms  or 
molecules.  Non-conserved  variables  include,  but  are  not  limited 
to,  order  parameters  associated  with  crystal  structure,  e.g.,  its 
symmetry  and/or  lattice  orientation.  The  present  work  deals 
only  with  non-conserved  variables,  i.e.,  order  parameters.  In  this 
context,  the  term  “phase”  denotes  a  certain  crystal  structure  or 
lattice  configuration.  In  regions  of  uniform  phase,  order  parameters 
take  on  discrete  values  of  zero  or  one.  In  interfacial  regions, 
order  parameters  enable  interpolation  between  pure  phases.  In 
addition  to  depending  on  usual  mechanical  and  thermal  state 
variables  (e.g.,  strain  and  temperature),  the  local  free  energy 
density  of  a  substance  generally  depends  on  local  value(s)  of 
order  parameter(s)  and  spatial  gradients  of  order  parameter(s). 
Such  a  prescription  enables  representation  of  surface  energies  of 
phase  boundaries  associated  with  order  parameter  gradients  in 
interfacial  regions. 
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Pioneering  treatments  of  thermodynamics  and  kinetics  of  het¬ 
erogeneous  material  systems  in  the  context  of  phase  field  models 
were  forwarded,  respectively,  by  Cahn  and  Hilliard  [  1  ]  and  Allen 
and  Cahn  [2].  A  fundamental  ansatz  [2]  often  prescribed  in  tra¬ 
ditional  phase  field  modeling  is  that  a  material  system  will  tend 
to  evolve  towards  a  state  of  minimum  free  energy,  subject  to 
boundary  constraints  imposed  on  the  system.  Concepts  for  mod¬ 
eling  multiphase  systems  were  advanced  by  Steinbach  et  al.  [3] 
and  Steinbach  and  Apel  [4].  Fried  and  Gurtin  [5]  and  Gurtin  [6] 
developed  order  parameter  theories  incorporating  configurational 
forces  or  micro-force  balances  in  the  context  of  geometrically  non¬ 
linear  continuum  mechanics  and  thermodynamics.  Review  articles 
are  available  that  describe  various  numerical  techniques  and  appli¬ 
cations  [7,8]. 

In  the  present  work,  phase  field  theory  is  used  to  describe  defor¬ 
mation  twinning,  also  known  as  mechanical  twinning,  defined  as 
twinning  induced  by  mechanical  stresses  [9-11].  Hereafter  defor¬ 
mation  twinning  is  simply  referred  to  as  “twinning”.  The  transfor¬ 
mation  strain  associated  with  twinning  is  a  simple  shear.  Across 
the  habit  plane,  the  orientation  of  the  Bravais  lattice  differs  by  a 
reflection  or  rotation  depending  on  the  kind  of  twin  under  consid¬ 
eration.  The  sheared  and  re-oriented  crystal  is  termed  the  “twin”, 
while  the  region  of  crystal  that  maintains  its  original  orientation  is 
termed  the  “parent”  or  the  “matrix”.  Twinning  often  takes  place  by 
coordinated  movement  of  partial  dislocations  (i.e.,  twinning  dislo¬ 
cations)  and/or  shuffles  of  some  but  not  all  atoms  comprising  the 
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twin.  A  finite  energy  can  be  associated  with  the  twin  boundary, 
typically  estimated  on  the  order  of  the  appropriate  stacking  fault 
energy  [12].  More  generally,  the  interfacial  energy  of  a  growing  or 
shrinking  twin  also  includes  elastic  and  core  energies  of  twinning 
dislocations  comprising  such  interfaces  [9,13,14]. 

Continuum  mechanics  models  for  deformation  twinning  have 
been  developed  in  the  context  of  crystal  plasticity  theory  [15,16], 
wherein  at  each  continuum  material  point,  the  volume/mass 
fraction  of  a  particular  twin  system  is  evolved  via  a  kinetic 
relation  usually  involving  a  resolved  shear  stress  criterion. 
Intended  goals  of  such  models  include  predictions  of  macroscopic 
stress-strain  behavior  and  crystallographic  texture.  These  models 
are  generally  unable  to  predict  detailed  twin  morphologies  at  the 
nanometer  scale.  Sharp  interface  models  have  been  developed 
to  address  kinetics  of  twin  growth  [17,18].  Such  models  are 
capable  of  predicting  twin  morphology,  but  application  of  these 
theories  requires  prescription  of  kinetic  laws  for  motion  of 
twin  boundary  interfaces,  and  their  numerical  solutions  require 
advanced  computational  methods  such  as  level  sets  [18]  for 
resolution  of  surfaces  of  discontinuity. 

Continuum  phase  field  approaches  have  been  suggested 
elsewhere  for  twin  growth  kinetics  [19,20]  invoking  the  Time- 
Dependent  Ginzburg-Landau  (TDGL)  approach  for  energy 
minimization.  A  different  approach  is  pursued  in  the  theory  and 
numerical  solution  techniques  developed  in  the  present  paper.  The 
general  theory  developed  here  addresses  the  following  physics: 
potential  activity  of  one  or  more  twin  systems,  large  deformations, 
nonlinear  elastic  behavior,  and  elastic  and  surface  anisotropy.  Or¬ 
der  parameter(s)  are  related  to  magnitude(s)  of  twinning  shear 
(i.e.,  transformation  strain  and  rotation)  at  a  material  point.  Inves¬ 
tigated  in  this  work  are  equilibrium  phenomena  in  the  null  tem¬ 
perature  (i.e.,  0  K)  and  quasi-static  limits  that  can  be  addressed 
via  energy  minimization  and  thermodynamic  stability  principles. 
Dissipation  and  time  scales  associated  with  growth  kinetics  are 
not  formally  addressed,  nor  are  acoustic  waves.  The  present  ap¬ 
proach  is  somewhat  analogous  to  lattice  statics  treatments  of  dis¬ 
location  mechanics  [21,22],  wherein  dislocation  slip  in  the  OK  limit 
is  addressed  without  consideration  of  atomic  vibrations  associated 
with  dissipated  heat  from  defect  motion.  Quasi-static  approaches 
have  similarly  been  used  to  study  twinning  in  the  context  of  empir¬ 
ical  pair  potentials  [23],  empirical  many-body  potentials  [24],  and 
density  functional  theory  (DFT)  [25].  Idesman  et  al.  [26]  developed 
a  quasi-static  continuum  approach  for  modeling  elastic-plastic 
materials  undergoing  martensitic  phase  transitions  and  twinning. 
Koslowski  et  al.  [27]  formulated  a  phase  field  theory  of  disloca¬ 
tion  mechanics  incorporating  a  non-convex  Peierls  potential  and 
energy  minimization  concepts.  The  advantage  of  the  present  ap¬ 
proach  over  TDGL  methods  is  that  material  parameters  associated 
with  time  scales  for  interfacial  motion  do  not  enter  the  model  and 
need  not  be  measured  experimentally.  This  is  an  important  consid¬ 
eration  for  modeling  of  deformation  twinning  since  the  propaga¬ 
tion  speed  of  twin  boundaries  can  be  difficult  to  measure,  and  could 
even  be  supersonic  if  the  driving  stress  is  sufficiently  large  [17,28], 
though  twin  propagation  speeds  in  the  subsonic  regime  have  also 
been  observed  [29].  Modeling  of  twin  growth  kinetics  at  realistic 
time  scales  in  the  former  case  would  seem  to  require  resolution  of 
stress  dynamics,  i.e.,  wave  propagation  and  possible  shock  wave 
phenomena. 

The  present  application  of  the  theory  considers  homogeneous 
twin  nucleation  in  an  otherwise  defect-free  single  crystal  subjected 
to  far-field  stress  [30-33].  For  given  far-field  boundary  conditions, 
the  minimum  stable  size  and  equilibrium  shape  of  a  twin  embryo 
are  dictated  by  competition  between  elastic  strain  energy  and  sur¬ 
face  energy  associated  with  the  twin-parent  boundary.  The  theory 
is  implemented  in  a  finite  element  code  that  seeks  minima  of  a  free 
energy  functional  that  depends  on  deformation  gradient  and  order 


parameter  fields.  When  the  surface  energy  is  idealized  as  isotropic, 
the  only  material  parameters  required  are  the  characteristic  twin¬ 
ning  shear  and  geometry  (known  from  the  crystal  structure),  the 
elastic  constants,  the  twin  boundary  energy,  and  a  characteristic 
length  associated  with  the  equilibrium  twin  boundary  thickness. 
Numerical  simulations  are  used  to  investigate  criteria  for  twin  nu¬ 
cleation  and  growth  depending  on  critical  resolved  shear  stress, 
twin  boundary  energy,  and  far-field  boundary  constraints.  Specif¬ 
ically  considered  in  this  paper  are  (1011)  {1012}  twins  in  magne¬ 
sium  (Mg)  single  crystals  [13,14,34,23]. 

The  remainder  of  this  paper  proceeds  as  follows.  Section  2 
presents  the  geometrically  nonlinear  phase  field  theory  for 
mechanical  twinning.  Section  3  presents  a  geometric  linearization 
of  the  theory,  wherein  deformations  are  assumed  small  as  in 
conventional  linear  elasticity.  The  model  described  in  Section  3 
is  still  nonlinear  in  the  sense  that  the  total  free  energy  of  the 
material  is  addressed  by  a  potential  non-quadratic  in  the  order 
parameter.  Section  4  presents  numerical  methods  used  to  obtain 
solutions  of  boundary  value  problems,  specifically  finite  element 
methods  with  free  energy  minimization  proceeding  via  a  conjugate 
gradient  algorithm.  Section  5  presents  model  predictions  for  twin 
nucleation  in  Mg  single  crystals.  Conclusions  follow  in  Section  6. 

Notation  of  continuum  mechanics  is  used  [5,6,10].  Real 
numbers  are  M.  Vectors  and  higher-order  tensors  are  written  in 
bold  font;  scalars  and  components  of  vectors  and  tensors  are 
written  in  italic  font.  When  indicial  notation  is  used,  summation 
proceeds  over  repeated  indices.  To  simplify  notation,  vectors  and 
tensors  are  referred  to  a  fixed  Cartesian  frame  of  reference,  with 
indices  in  the  subscript  position.  The  scalar  product  of  vectors 
a  and  b  is  a  •  b  =  =  Oibi  -h  02^2  +  03^3  in  a  three- 

dimensional  vector  space.  The  outer  product  is  (a  0  bj^e  = 
Juxtaposition  implies  summation  over  one  set  of  adjacent  indices, 
e.g.,  (AB)ab  =  AacBcb-  The  colon  denotes  summation  over  two 
sets  of  indices;  e.g.,  A  :  B  =  AabBab  and  (C  :  E)ab  =  CabcdEcd-  The 
transpose  of  a  matrix  is  indicated  by  a  T  superscript,  e.g.,  =  Aba- 

The  inverse  operation  is  denoted  by  a  —1  superscript,  and  the 
inverse-transpose  operation  is  denoted  by  a  —T  superscript,  e.g., 

2.  Geometrically  nonlinear  theory 

In  what  follows  in  Section  2,  a  phase  field  theory  is  developed 
for  an  elastic  body  with  a  single  twin  system  and  describable 
by  a  single  order  parameter.  Extension  of  the  theory  to  elastic 
bodies  with  multiple  twin  systems  and  multiple  order  parameters 
is  considered  in  Appendix  B. 


2.1.  Order  parameter 


Let  ^2  c  be  a  reference  configuration  of  a  body  and  X  g 
^2  be  a  material  point.  Existence  of  an  order  parameter  function 
77  :  i?  X  (0,  T)  ^  [0,  1],  where  (0,  T)  c  M  is  a  time  interval, 
is  assumed.  The  order  parameter  function  distinguishes  between 
two  distinct  phases:  (i)  the  original  crystal  (the  parent)  and  (ii)  the 
twin.  Interfaces  between  phases  represent  twin  boundaries.  Order 
parameter  p  exhibits  the  following  values: 


^(X,  •) 


0  ifXG  parent, 

(0,1)  ifXG  twin  boundary, 
1  ifXG  twin. 


(1) 


According  to  diffuse  interface  theory  [1,2],  p  is  commonly 
presumed  at  least  continuous  with  respect  to  X. 


2.2.  Kinematics 

A  motion  of  ^2  on  (0,  T)  c  M  is  given  by  a  map  x  :  ^2  x  (0,  T)  ^ 
K^.  Spatial  coordinates  x  and  reference  coordinates  X  e  of  a 
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material  particle  are  thus  related  by 

x=X(X,  t).  (2) 

The  deformation  gradient  is 

F  =  Vx,  FaA  =  VAXa,  (3) 

with  =  d/dXA  the  material  gradient.  In  regions  of  the  body 
where  x  is  at  least  with  respect  to  X,  associated  compatibility 
conditions  are  V  x  F  =  V  x  (Vx)  =  0  where  x  is  the  vector 
cross  product.  Associated  with  twinning  kinematics  is  the  simple 
shear [10] 

F  =  1  +  yoS  0  m  (with  s  •  m  =  0  and  s  •  s  =  m  •  m  =  1).  (4) 

The  unit  tensor  is  1.  The  unit  normal  to  the  surface  of  composition 
(i.e.,  the  habit  plane)  in  the  reference  configuration  is  m.  The 
magnitude  of  twinning  shear  and  shear  direction  are  yo  and  s. 
respectively;  these  are  both  constants  in  a  homogeneous  crystal. 
Let  F+  and  F“  denote  deformation  gradients  in  twin  and  parent, 
respectively,  both  in  their  null  strain  energy  reference  states.  The 
difference  between  these  deformation  gradients  is  the  rank-one 
matrix 

F+  —  F“  =  F  —  1  =  yoS  0  in.  (5) 

An  interpolation  procedure  [20]  is  used  to  characterize  the 
twinning  shear  in  the  interfacial  regions.  Deformation  gradient  (3) 
is  decomposed  multiplicatively  as 

F  =  F^F'', 
where 

F^(X,  t)  =  F(F'^)“^  =  elastic  deformation, 

F'^[r](X,  t)]  =  stress-free  twinning  shear. 

Superscripts  E  and  r]  are  descriptive  labels  and  not  numerical  ex¬ 
ponents.  Specifically,  twinning  shear  is  interpolated  in  interfacial 
regions  as  follows: 

F^  =  1  -h  [(p(r])]  yo  s  (g)  m,  (8) 

where  tp  :  [0,  1]  ^  [0,  1]  is  a  monotonically  increasing  function 
obeying  ^(0)  =  0  and  (p{\)  =  1.  A  representative  function  also 
satisfying  the  condition  of  vanishing  derivative  at  the  endpoints 
[^'(0)  =  (p\\)  =  0]  that  will  be  used  later  is  a  “2-3-4 

polynomial”  [35]. 

(p{r])  =  arf'  +  2(2  —  -h  (o'  —  (9) 

where  a  is  a  scalar  constant  within  the  limits  0  <  a  <  6.  It 
follows  that  in  the  twin,  F^(l)  =  1  -h  yoS  (g)  m  =  F,  and  in 
the  parent,  F^(0)  =  1.  Twinning  preserves  the  volume  and  mass 
density  of  a  material  element  of  the  solid:  detF"^  =  det(l  -h 
^yoS  (g)  m)  =  1  -h  ^yoS  •  m  =  1.  Fig.  1  shows  an  elastic 
body  undergoing  twin  nucleation.  The  intermediate  configuration, 
twinned  but  elastically  unloaded,  is  labeled  “fictitious”  as  a  result 
of  incompatibility  conditions  V  x  F^  /  0  [36].  While  the  present 
theory  is  fully  three-dimensional,  a  two-dimensional  coordinate 
system  is  inscribed  on  the  reference  body  for  later  reference  in 
Section  5.  The  theory  presented  here  does  not  address  plastic 
slip,  i.e.,  glide  of  full  and/or  partial  dislocation  lines  and  loops  not 
associated  with  twinning.  Additional  model  features  are  required 
to  simultaneously  address  plastic  slip  and  twinning  [16,36]. 

2.3.  Free  energy  functional 


(6) 

(7) 


\  / 
\  / 


Fig.  1.  Phases  and  kinematics  of  elasticity  and  twinning. 


where  W  :  i?  x  x  [0,  1]  ^  R  is  the  elastic  strain  energy 
density  (generally  nonzero  in  parent,  twin,  and  interfacial  regions), 
and  where/  :  i?  x  [0.1]  X  R^  ^  R  accounts  for  interfacial  energy 
in  twin  boundary  regions.  R+^^  denotes  the  set  of  3  x  3  matrices 
with  positive  determinant.  Henceforth,  the  following  functional 
form  of  the  strain  energy  per  unit  reference  volume  is  used: 

W(X,  F,  r])  =  W[E^(F,  rj),  rj], 

E^  =  l(C"-l)  =  l(F"V-l),  (11) 

with  the  elastic  deformation  tensor  and  the  elastic  strain 
tensor.  The  elastic  strain  can  be  expressed  in  terms  of  F  and  r]  via 
use  of  (7)  and  (8).  For  any  value  of  order  parameter  strain  energy 
density  vanishes  at  null  elastic  strain: 

W(0,  •)  =  0.  (12) 

A  quadratic  form  for  the  strain  energy  density  is  assumed. 
Higher-order  elastic  coefficients  [16,20,37]  can  be  incorporated 
by  straightforward  extension.  Dependence  of  W(E^,  ri)  on  r] 
manifests  explicitly  only  via  anisotropic  elastic  coefficients.  Strain 
energy  density  and  second-order  moduli  are  written 


W  =  :  C(r])  :  E"^ 


C(ri)  = 


d^w 

dE^dE^ 


E^=0 


(13) 


As  usual,  indices  Cabcd  =  Qadc  =  Qdab  ^nd  the  Laue 
group  of  the  crystal  dictates  any  other  symmetries  of  C(0) 
in  the  crystallographic  frame  [37].  For  a  compound  twin  in  a 
centrosymmetric  structure,  re-orientation  matrix  Q  associated 
with  twinning  is  [  1 1  ] 


(i=  2m  (g)  m  —  1.  (14) 

Elastic  coefficients  of  the  fully  twinned  crystal  are  related  to  those 
of  the  parent  by 


The  total  free  energy  functional  for  a  body  undergoing  twinning 
deformation  is  written  as 

I  W(X,E,r])dQ  [  f  (X,  r] ,  V r])  d^2 ,  (10) 

J  £2  J  £2 


Qbcd(^)  =  QaeQbfQcgQdhQfgh(O).  (15) 

Elastic  coefficients  in  interfacial  regions  are  interpolated  the  same 
way  as  the  twinning  shear: 

C(rf)  =  C(0)  + [CO) -C(0)]  (p(rf), 


(16) 
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where  (p(r])  obeys  (9).  Different  interpolators  could  be  used  for 
twinning  shear  and  elastic  coefficients;  for  simplicity,  the  same 
function  is  prescribed  here  for  both.  In  the  isotropic  elastic 
approximation,  letting  X  denote  Lame’s  constant  and  /x  the  shear 
modulus,  the  elasticity  tensor 

CaBCD  =  ^^AB^CD  +  f^i^AC^BD  +  ^AD^Bc),  (^7) 

SO  that  C(0)  =  C(l)  and  W(E^,  ^  W(E^)  does  not  explicitly 

depend  on  r]. 

The  local  interfacial  energy  per  unit  reference  volume  follows 
from  the  Cahn-Hilliard  formalism  [1]: 

fill,  =foir])  -\-ic  :  (Vt]  <S>  V/y),  (18) 

with  K  a  symmetric  second-order  tensor  that  may  generally 
depend  on  rj,  but  is  assumed  here  for  simplicity  to  have  constant 
components.  When  interfacial  energy  is  isotropic,  k  =  k  \  and 

f{ri,Vvi)={o{ri)+K\V^t  (19) 

Prescribed  for/o  is  a  standard  “double-well”  potential  [5,7,20,35]: 

foirii)=Arf{\-rif,  (20) 

with  A  >  0.  As  shown  in  Appendix  C,  in  the  isotropic 

approximation  A  and  k  are  related  to  equilibrium  energy  per  unit 
area  r  and  thickness  I  of  an  unstressed  interface  via 

K  =  3r//4,  A  =  nr/x.  (21) 

The  total  free  energy  functional  ^  of  (10)  becomes,  using  (13)  and 

(18). 

=  l  f  :  Cin)  :  E^di2 

+  [  [At]^(\-t]y  +  K:iyt]®yt])]dQ.  (22) 

JQ 

For  isotropic  elastic  and  interfacial  energies,  with  trA  =  Aaa  the 
trace  of  a  second-order  tensor, 

ri)=  f  [(X/2)(trE")2  +  :  E"]di2 

+  I  [At]^(\-r]y  +  K\yt]\^]di2.  (23) 

J  C2 

2.4  Equilibrium  conditions 


order  parameter.  Assuming  -smoothness  off  and  taking  the  first 
variation  of  the  interfacial  energy, 

8  f  d^2  =  f' 8r]  d^2 2ic  :  (Vr]  (g)  V87])  d^2,  (26) 

J  r2  J  J  r2 

where  8r]  denotes  an  admissible  variation  of  order  parameter 
rj.  Similarly,  for  W  a  -function  of  all  its  arguments,  the  first 
variation  of  the  strain  energy  is 


r  f  dW  f  dW 

8  Wd^=  - 8r]d^-\-  -  :  V8x  dQ,  (27) 

JQ  Jc2  ^11  JQ 

with  8x  denoting  an  admissible  variation  of  deformation  map 
X.  Combining  (25)-(27)  yields  the  weak  form  of  the  equilibrium 
equations: 

“h  — — ^  8t]  dif2  “h  2ic  :  (V/y  (g)  dif2 

-h  /*  —  :  V5xd^2  -  /*  t-8xdS-f  lj8r]dS  =  0  (28) 

JQ  JdQ  JdQ 

for  any  admissible  variations  8x  and  8ri. 

Strong  forms  of  equilibrium  equations  require  C^-smoothness 
of  both  X  and  q.  Application  of  the  divergence  theorem  to  (26)  leads 
to 

8  fd^  =  f'8r]d^2  -  2ic  :  [V(V^)]  8r]  dQ 

J  if?  J  if?  J  if? 


-h  /  2ic  :  [(Vq)  (g)  n]  8q  dS,  (29) 

JdQ 

with  n  the  unit  outward  normal  to  9^2.  The  first  variation  of  the 
strain  energy  can  be  converted  analogously  into 

S  f  Wd^2  =  f  —St^d^2  -  f 

JQ  JQ  ^9  JQ 


dw 


■  8x  d^2 


(30) 


Local  strong  forms  of  equilibrium  equations  (Euler-Lagrange 
equations)  follow  directly  from  (29)  and  (30): 


V  • 


dw 


VP  =  0, 


dw 


dVAXa 


^aPqa  =  0, 


(31) 


A  stable  configuration  of  a  body  undergoing  twinning  deforma¬ 
tion  corresponds  to  a  minimizer  of  total  free  energy  functional  (10), 
given  certain  boundary  conditions.  The  problem  of  finding  such 
stable  configuration  can  be  expressed  as 

mm^(X,q).  (24) 

X,r] 

Substantial  literature  [38,39]  has  been  dedicated  to  assessment 
of  existence  of  minimizers  such  as  in  (24).  In  general,  such 
minimizers  need  not  exist  for  arbitrary  energy  density  functions 
/  and  W.  Moreover,  existence  requirements  impose  constraints 
on  functional  forms  of  /  and  W.  In  Appendix  A,  existence  of 
minimizers  of  (24)  is  explored  under  assumptions  of  convexity  or 
polyconvexity  of  strain  energy  density  function  W  and  interfacial 
energy  density  function/  [40,39]. 

The  following  variational  equation  is  posited  that  will  suggest, 
upon  application  of  Hamilton’s  principle,  both  weak  and  strong 
forms  of  static  equilibrium  equations  and  boundary  conditions: 

8^-1  t  •  9x  d5  -  /  h8qdS  =  0,  (25) 

JdQ  JdQ 

where  t  is  a  mechanical  traction  vector  per  unit  reference  area,  9^2 
is  the  boundary  of  ^2,  d5  is  a  surface  element  of  9^2,  and  h  is  a  scalar 
conjugate  force  (energy  per  unit  reference  area)  to  variations  of  the 


(32) 


,  9W 

f'  -  2k  :  [V(V/7)]  +  ^  =0, 

9ri  F 

9VV 

fo  —  '2'KABy Bf]  +  —  —  0, 

dtl  F 

where  P  is  the  first  Piola-Kirchhoff  stress  tensor.  Corresponding 
boundary  conditions  are 

t  =  Pn,  ta  =  PaAnA,  (33) 

h  =  2k  :  (Vri  (gin),  h  =  iKAsns^ a^I ■  (34) 

The  first  Piola-Kirchhoff  stress  also  obeys 


dW 


F'=X'(F’>) 


dW  _  9E^  _  9F^ 

The  elastic  second  Piola-Kirchhoff  stress  is 
9W 


n\-T 


C  :  E'= 


(35) 


(36) 


The  partial  derivative  of  W[E^(F,  r]),  r]]  with  F  fixed  is  computed 
as 


9W 

dw 

9W 
-j-  — — 

9E^ 

dq 

F  ^9 

fe  9E£ 

^  '  9'? 

(37) 
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From  (13)  and  (16), 
dW 
drj 


1  P  dc  p 

=  -E^  :  —  : 


dr] 

l^E^  :  [C(l)  -  C(0)]  :  E^ 
2  drj 


(38) 


From  (6)-(8)  and  (36), 


aw 


dE^ 

drj 


-{X:  [d(s( 


dtp 

drj 


rniXF^r’llKo^ 

drj 


(39) 


It  can  be  verified  that  r  is  a  resolved  shear  stress  on  the  habit  plane 
in  the  direction  of  twinning  shear.  Combining  (37)-(39),  phase 
equilibrium  condition  (32)  can  be  rewritten 

1 


f'  -  2ic  :  [V(Vr])]  =  -E"  :  [C(0) 


■C(l)]:E"  +  rKoj^.(40) 

In  the  isotropic  approximation,  choosing  cp  from  (9)  and  using /o 
from  (20),  equilibrium  condition  (40)  becomes 

x[arf  +  3(2  —  a)rf‘  +  2{a  —  3)?]^] 

= —[At](\-3t]  +  2r]^)-KV^t]].  (41) 

Ko 

Both  sides  of  (41)  vanish  identically  in  regions  of  uniform  phase 
where  rj  =  0  or  rj  =  1.  For  purposes  of  comparison,  a  kinetic 
equation  corresponding  to  (40)  in  the  context  of  TDGL  theory  is 
listed  in  Appendix  D. 

3.  Geometric  linearization 

In  what  follows  in  Section  3,  the  theory  of  Section  2  is 
linearized  for  small  deformations.  The  linearized  theory  proves 
useful  for  comparison  of  numerical  results  with  analytical  studies 
of  energetics  of  twin  nucleation  [30-33]  in  the  context  of  Eshelby’s 
model  of  elastic  inclusions  and  inhomogeneities  [41 ,42  ].  Eq.  ( 1 )  still 
applies. 

3A.  Kinematics 

In  linear  elasticity,  the  usual  kinematic  field  variables  are 
displacement  u(X,  t)  =  x(X,  t)  —X  and  its  gradient 

p  =  Vu,  Pab  =  ^bUa.  (42) 

In  contrast  to  the  geometrically  nonlinear  theory  (Fig.  1),  there 
is  no  explicit  distinction  among  configurations  of  a  deformable 
body.  Compatibility  conditions  are  V  x  ^  =  V  x  (Vu)  =  0. 
Decomposition  (6)  is  replaced  with 

P  =  K  +  (43) 

where  is  the  elastic  distortion  and  is  the  distortion  associated 
with  twinning  shear: 

=  [(p(r])]yo  s  (g)  m.  (44) 

The  symmetric  elastic  strain  tensor  is 


1 


-{Vu  +  (Vu)^  -  yo[(p{r])][s  (g)  m  +  m  (g)  s]}.  (45) 


3.2.  Free  energy 

The  total  free  energy  functional  for  a  body  of  reference  volume 
F2  is  written  as 

r])=  W(Vu,  r])d^2  +  /  f(r],  Vr])d^2, 

Jr?  Jr? 


(46) 


where  W  is  the  elastic  strain  energy  density  and  /  accounts  for 
interfacial  energy.  Strain  energy  density  and  second-order  elastic 
moduli  are 


w  =  W[e^(Vu,  rf),  /?]  =  :  C(j))  : 


€(/?)  = 


d^W 


de^de^ 


(47) 


e^=0 


Eqs.  (14)-(21)  apply  in  the  linearized  case.  The  total  free  energy 
functional  ^  of  (46)  becomes,  using  (18)  and  (47), 

^ (u,  rj)  =  -  f  :  C(r])  :  di? 

2  Jr? 


/i 


+  I  [Ari^ -  ri)^  +  K  :  (V ri  ®  y  ri)]  dQ . 


(48) 


For  isotropic  elastic  and  interfacial  energies,  this  reduces  to 

>I'(u,t])=  /  [(X/2)(trVu)^  +  /r(Vu-Ko?'S(8)in)symm 

Jr? 

:  (Vu  -  Ko?' s  (8>  in)symm]  df2 


+  I  [At]^(\-r]Y+K\yr]\^]dQ, 

Jr? 


(49) 


where  (•)symm  denotes  the  symmetric  part  of  a  second-order 
tensor,  e.g.,  2Asymm  =  A  -h  A^. 

3.3.  Equilibrium  conditions 

The  following  variational  equation  is  posited  that  will  suggest, 
upon  application  of  Hamilton’s  principle,  local  or  strong  forms  of 
static  equilibrium  equations  and  boundary  conditions: 


8^  -  t -^udS  -  / 

Jdr?  Jdr? 

where  t  is  a  mechanical  traction  vector  per  unit  area,  dS  is  a  surface 
element  of  9^2,  and  h  is  a  scalar  conjugate  force  to  variations  of  the 
order  parameter.  Taking  the  first  variation  of  the  interfacial  energy 
and  applying  the  divergence  theorem  gives  (29).  The  first  variation 
of  the  strain  energy  with  Vu  and  rj  independent  is 

o  /*  ,  /*  9^0  ,  /*  r_  9W  ^ 

8  Wd^2  =  /  - 8r]dQ  - 

Jr?  Jr?  Jr? 


h8r]dS  =  0, 


(50) 


V  • 


9Vu 


■8ndQ 


+ 


/  [" 


dw ' 

Wix 


■  8u  dS. 


(51) 


It  follows  from  (29)  and  (51)  that  in  ^2,  the  Euler-Lagrange 
equations  are 


V  • 


dw 

d^ 


=  V  a  =  0, 


(52) 


^  0, 


,  dW 

f'  -  2fc  :  [V(Vr])]  +  — 
drj 

where  a  is  the  symmetric  stress  tensor.  Corresponding  boundary 
conditions  are 

t  =  an,  h  =  2ic  :  (Vrj  (g)  n). 

The  stress  also  obeys,  from  (47), 


a  = 


dw 

dVn 


dW  de^ 
de^  ’  9Vu 


(53) 

9W  F 

—  =  C  : 

de^ 

(54) 

The  partial  derivative  ofW[e^  (Vu,  rj),  rj]  with  Vu  fixed  is 

dW 
drj 


dW 

dw 
— — 

de^ 

Vu 

dr] 

(55) 
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From  (16)  and  (47), 


dW 

drj 


-e 

2 


dc 

drj 


:  [C(l)  -  C(0)]  : 

2  drj 


From(42)-(45)and(54), 


dW 

de^ 


de^ 

dr] 


dtp  dtp 

=  -{or  :  (s  (8)  m)}yo^  = 

drj  drj 


(56) 


(57) 


where  r  is  a  resolved  shear  stress  acting  on  the  habit  plane  in  the 
direction  of  twinning  shear.  Combining  (55)-(57),  the  second  of 
equilibrium  conditions  (52)  can  be  rewritten 


/o'  -  2k  :  [V(V;,)]  =  :  [C(0)  -  C(l)]  :  +  ryoj  (58) 

In  the  isotropic  approximation,  choosing  tp  from  (9)  and  fo 
from  (20),  (58)  reduces  to  (41),  but  with  r  defined  in  (57)  rather 
than  (39).  When  elastic  strains  are  small,  X"  ~  or,  ~  1,  and  the 
difference  between  r  in  (39)  and  (57)  is  on  the  order  of  twinning 
shear  yo. 


4.  Numerical  methods 


The  finite  element  method  is  used  to  seek  solutions  for 
equilibrium  or  local  minimum  energy  states  of  a  body  subjected  to 
boundary  conditions.  At  each  reference  point  X,  primary  solution 
variables  are  displacement  u  and  order  parameter  rj.  If  these 
solution  variables  and  requisite  material  properties  are  known, 
then  all  elastic  field  variables  and  interfacial  quantities  can  be 
computed  via  the  appropriate  mathematical  operations  given  in 
Sections  2  and  3. 


4. 2.  Finite  element  discretization 


The  body  of  reference  volume  Q  is  discretized  into  a  number 
of  standard  finite  elements  with  shape  functions  Ni(X).  Let 
U/(t)  and  r]i(t)  denote  instantaneous  values  of  displacement  and 
order  parameter,  respectively,  at  node  i.  Displacement  and  order 
parameter  fields  are  represented  in  a  finite  element  context  as, 
respectively, 

u(X,  t)  =  Ui(t)Ni(X),  r](X,  t)  =  r]i(t)Ni(X),  (59) 

where  summation  proceeds  over  all  nodes  i;  however,  only  those 
nodes  supporting  the  element(s)  containing  or  bounding  point  X 
have  Ni(X)  /  0.  Material  gradients  of  (59)  follow  as 

F  =  1  +  Vu  =  1  +  Uj  (g)  VNj,  Vr]  =  rjiVNi.  (60) 


4.2.  Energy  minimization 

The  finite  element  method  is  employed  to  seek  solutions  of 
weak  forms  of  equilibrium  conditions  (  28).  Strong  forms  derived  in 
Sections  2.4  and  3.3  are  not  needed  by  the  numerical  algorithms. 
Addressed  in  what  follows  are  the  following  kinds  of  boundary 
conditions: 

9i7  =  ^  dF2p 

d^M  —  dF2M,D  U  dQM,N’)  0  —  dQM,N 

u(X,  t)  prescribed  on  9^2m,d,  t(X,  t)  =  0  on  9^2m,n 

mechanical  conditions  (61) 

9i2p  =  d^2p  Q  U  d^2pj\j,  0  =  d^2p  q  H  d^2p  j\f 
r](X,  t)  prescribed  on  dQp^o,  h(X,  f)  =  0  on  d^2p^pi 

phase  field  conditions 

Dirichlet  conditions  for  displacement  and  the  order  parameter 
are  applied  on  9^2m,d  and  9^2p,d.  respectively.  Neumann  condi¬ 
tions  corresponding  to  free  surfaces  are  applied  on  9^2m,n  and 


9^2p,Af,  respectively.  The  present  treatment  can  easily  be  extended 
to  address  arbitrary  Neumann  conditions  if  terms  accounting  for 
external  work  are  incorporated  in  the  energy  functional  whose  sta¬ 
tionary  points  are  sought. 

First  consider  the  geometrically  nonlinear  theory  of  Section  2. 
The  free  energy  functional  is  (22)  in  the  general  anisotropic  case. 
Nodal  equilibrium  conditions  are  obtained  by  substituting  (59) 
and  (60)  into  (22)  and  differentiating  with  respect  to  nodal  degrees 
of  freedom: 


_  f 

dUi  j Q 

/ 

Jr2 


dW  9F 

.  -  :  —  d^2 

dUi  Jq  9F  dUi 


/ 

Jr2 


PVNi 


[F^(C  :  E^)(F^)-^]VNj  d^2, 


(62) 


0  = 


a 

drji 


f  —Ni  d^2+  f  -2—  ■  VNi  +  /  — N; 
Ja  dri  Ja  9V??  Jn  dr]  ' 


9/ 


=  / 


+ 


L\'2^‘ 


/ 

Jr2 


dw 


{2Ar](l  -3r]-F  2r]^)}Ni  dQ  -h  /  (2/rV/y)  •  VN,-  d^2 

J  r2 

:  [C(l)  -  C(0)]  :  -  T/ol  d^2.  (63) 

I  dr] 


In  general,  solutions  of  (62)  and  (63),  corresponding  to  stationary 
points  8^(\Xi,rfi)  =  0,  can  be  associated  with  local  energy 

minima,  maxima,  or  saddle  points.  Energy  functional  ^(u/,  r]i)  is 
not  necessarily  convex  in  its  arguments,  and  may  exhibit  multiple 
local  minima,  for  example  [20].  A  conjugate  gradient  algorithm  is 
used  to  seek  minimum  energy  states  corresponding  to  (62)  and 
(63).  Initial  conditions  are  also  prescribed  as  part  of  the  solution 
procedure.  For  example,  in  Section  5  an  initial  twin  nucleus  {r]  =  \) 
is  placed  within  a  larger  domain  wherein  t]  =  0  initially.  Such  a 
system  will  not  be  in  mechanical  or  phase  equilibrium  at  the  initial 
time. 

Now  consider  the  linearized  model  of  Section  3,  with  energy 
functional  (48).  Global  discretized  equilibrium  conditions  are 


d^  f 
^  ~  dUi  ~  A 
=  /(C: 

JQ 

J 


dW  _  9Vu 
^  9Vu  ■  dUi 
e^)VNid^, 


dQ 


=/ 

Jr2 


aVNi  di? 


(64) 


[2Ar]('l  -3r]  +  2t]^)]Ni  AQ 


Q 

dr]i  jQ 

+  /  (2a:V?j)  •  VNj  di2 
Jn 


+ 


L\P 


d(p 


[C(l)-C(0)]:e"-rKohr'^>d^-  (6^) 

J  dr] 

For  the  isotropic  case  (49)  with  interpolator  tp  from  (9),  these 
reduce  to 


0  = 


d^ 

9Ui 


/ 

JQ 


P(V  •  u)l  +  2/i[Vu  -  Yoiarf  +  2(2  - 


-j- (of  3)/^  ]s  0  in]synini}VN,- d^2, 
d4> 
dr]i 


(66) 


/ 

Jr2 


/ 


{2At]i-[  -3r]  +  2r]^)}Ni  dQ+  [2KVr]}  •  VN;  di2 

r2  J  r2 


+  /  {mKo[[o;j)^  +  2(2  -  a)r]^  +  (o;  -  3)r]'']yo 


-  2Vu  :  (s  <S>  rnjsymm] 
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X  [2ari  -h  6(2  —  a)ri  -h  4(q'  —  3)r]  ]}Ni  dQ. 


(67) 
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Coupling  among  r],  Vrj,  and  Vu  depicts  interaction  of  the  order 
parameter  and  its  gradient  with  elastic  fields. 

5.  Application:  twin  nucleation  in  magnesium  single  crystals 

The  phase  field  theory  and  numerical  methods  discussed  in 
Sections  2-4  are  applied  to  the  problem  of  homogeneous  twin 
nucleation. 

5.2.  Background:  homogeneous  twin  nucleation 

Homogeneous  twin  nucleation  is  defined  as  the  nucleation  of 
a  twin  embryo  within  an  otherwise  perfect  single  crystal  [11]. 
Analytical  models  based  on  free  energy  variation  concepts  in  the 
context  of  phase  transformations  [42]  have  been  applied  to  de¬ 
scribe  twin  nucleation  [11,30-33,43].  Such  analytical  approaches 
consider  nucleation  of  a  twin  embryo  of  idealized  geometry  -  an 
elliptical  cylinder  in  two  dimensions  or  an  ellipsoid  in  three  dimen¬ 
sions  -  embedded  in  an  infinite  medium,  with  a  perfectly  bonded 
(i.e.,  coherent)  sharp  interface  separating  inclusion  from  surround¬ 
ing  medium  (i.e.,  the  matrix).  The  solution  technique  involves  si¬ 
multaneous  solution  of  two  equilibrium  equations  associated  with 
stationary  points  of  the  total  (Gibbs)  free  energy  change  associated 
with  twinning.  These  two  equations  yield  the  critical  size  and  as¬ 
pect  ratio  of  the  inclusion,  as  outlined  in  Appendix  E.  Apparently, 
exact  solutions  are  available  only  for  approximations  of  linear  elas¬ 
ticity,  isotropic  surface  energy  associated  with  the  twin  boundary, 
and  traction  boundary  conditions  at  infinity.  The  solution  provides 
the  critical  size  and  shape  of  an  inclusion  for  a  given  set  of  material 
properties  -  elastic  constants,  twin  boundary  surface  energy,  and 
twinning  transformation  shear  -  and  far-field  applied  stress.  At  the 
critical  aspect  ratio,  the  critical  size  corresponds  to  unstable  equi¬ 
librium  [31,32],  i.e.,  a  saddle  point  on  the  free  energy  surface.  At  a 
given  far-field  applied  stress,  a  twin  nucleus  smaller  than  the  criti¬ 
cal  size  will  tend  to  shrink  and  disappear,  while  one  larger  than  the 
critical  size  will  tend  to  grow  in  an  unstable  manner.  The  larger  the 
applied  stress  component  resolved  on  the  habit  plane  in  the  twin¬ 
ning  direction,  the  smaller  the  critical  size,  meaning  that  nucle¬ 
ation  of  a  small  twin  becomes  more  energetically  favorable  as  the 
resolved  shear  stress  increases.  In  a  real  heterogeneous  material  at 
finite  temperature,  nuclei  of  various  sizes  and  shapes  emerge  as  a 
result  of  local  statistical  fluctuations  [11,31,32].  According  to  the 
theory,  those  nuclei  larger  than  the  critical  size  would  form  twins 
that  grow  until  interactions  with  other  defects,  grain  boundaries, 
or  external  surfaces  occur. 

In  the  present  paper,  the  phase  field  approach  is  used  to 
model  the  problem  of  homogeneous  twin  nucleation.  Numerical 
results  obtained  from  the  phase  field  model,  under  assumptions 
of  geometric  linearity  (Section  3)  and  plane  strain  boundary 
conditions,  are  compared  with  the  analytical  solution  [30]. 
Such  an  exercise  provides  validation  for  the  phase  field  model 
and  numerical  methods  advanced  here.  Additional  numerical 
simulations  consider  effects  of  anisotropic  twin  boundary  (surface) 
energy,  variable  equilibrium  thickness  of  the  twin  boundary 
region,  various  habit  plane  directions  (i.e.,  lattice  orientations), 
and  various  boundary  conditions  (Dirichlet  vs.  Neumann,  shear  vs. 
tension,  and  domain  shapes).  Such  additional  factors  cannot  all  be 
addressed  via  known  analytical  elasticity  solution  techniques. 

5.2.  Material 

Pure  magnesium  single  crystals  are  studied.  Properties  are 
listed  in  Table  1  with  supporting  references.  Magnesium  exhibits 
a  hexagonal  crystal  structure  and  is  centrosymmetric,  with 
c/a  =  1.6235.  Jhe  twinning  system  of  consideration  is  the 

primary  one:  {1011){1012}  with  twinning  shear  yo  =  (3  — 
c^/a^)/(3^/^c/a)  =  0.1295  [11,34].  All  five  independent  second- 
order  elastic  constants  [44]  are  listed  in  Table  1,  though  in 


Table  1 

Properties  for  pure  Mg  single  crystals. 


Parameter 

Value 

Definition 

Reference 

c 

5.200  A 

lattice  parameters  (0  K) 

[13] 

a 

3.203  A 

Vo 

0.1295 

shear  for  (1011){1012}  twin 

[11] 

Cu 

63.5  GPa 

second-order  elastic  constants: 

[44] 

C33 

66.5  GPa 

extrapolated  from  4.2  to  0  K 

Cu 

25.9  GPa 

Cl3 

21.7  GPa 

C44 

18.4  GPa 

X 

24.0  GPa 

Lame  constant  (Voigt  average) 

/X 

19.4  GPa 

shear  modulus  (Voigt  average) 

K 

36.7  GPa 

bulk  modulus  (Voigt  average) 

V 

0.276 

Poisson’s  ratio  (Voigt  average) 

r 

117mJ/m2 

twin  boundary  energy 

[13] 

i 

1.0  nm 

equilibrium  boundary  thickness 

[20,45] 

K 

0.0878  nj/m 

gradient  energy  parameter 

Eq.  (21) 

A 

1.404  GPa 

double-well  energy  parameter 

Eq.  (21) 

calculations  that  follow,  Voigt  averages  [46]  are  used  for  isotropic 
elastic  constants.  Magnesium  single  crystals  are  not  strongly 
anisotropic  elastically  [47]  [Cn  ~  C33,  Cu  ~  C13,  and  C44  ~  (Cn  — 
Ci2)/2],  so  the  isotropic  elastic  approximation  appears  reasonable. 
Reuss  averages  [46]  for  shear  and  bulk  moduli  are  19.3  GPa  and 
36.7  GPa,  respectively,  nearly  identical  to  Voigt  averages  in  Table  1. 

Twin  boundary  surface  energy  is  obtained  from  DFT  calcula¬ 
tions  [13].  In  most  phase  field  simulations  discussed  later,  surface 
energy  is  assumed  isotropic,  following  usual  theoretical  stud¬ 
ies  [30-33],  but  in  one  case  anisotropy  of  the  surface  energy  is  con¬ 
sidered.  In  a  material  coordinate  system  with  axes  aligned  parallel 
to  twinning  direction  and  habit  plane  normal  (Xi  ||  s  and  X2  ||  m), 
the  gradient  coefficient  entering  (48)  is  written 


In  the  isotropic  case,  ku  =  ^22  =  k.  The  anisotropic  case 
considered  later  is  /cii/2  =  2k22  =  k,  which  favors  boundaries 
extended  parallel  to  the  habit  plane.  Recall  from  (C.16)  that  surface 
energy  F  oc  The  rationale  for  anisotropic  surface  energy  is 
that  twinning  dislocations  at  a  moving  portion  of  the  boundary 
(i.e.,  an  incoherent  interface  in  the  terminology  of  Kosevich  and 
Boiko  [9])  would  contribute  core  and  elastic  energies  to  the  total 
surface  energy  of  the  interface.  In  contrast,  the  fully  formed 
(i.e.,  coherent)  twin  boundary  surface  would  have  less  energy 
than  such  a  moving  portion  because  it  does  not  contain  energy  of 
dislocations,  just  stacking  fault  energy  associated  with  reflection 
of  the  lattice  across  the  twin  boundary.  The  anisotropic  case 
considered  here  will  cause  the  twin  to  elongate  in  the  direction 
of  s  and  shorten  in  the  direction  of  m,  in  order  to  decrease  the 
contribution  of  to  the  energy  in  (48). 

The  equilibrium  thickness  I  over  which  atoms  deviate  from 
their  jdeal  positions  is  taken  as  1  nm,  corresponding  to  about 
five  {1012}  planes.  This  value  follows  from  theoretical  studies  of 
perturbed  atomic  coordinates  in  twinned  hexagonal  close-packed 
structures  [45].  The  same  characteristic  thickness  value  (1  nm)  has 
been  used  in  phase  field  models  of  other  materials  [20].  In  one  set 
of  calculations  that  follow,  a  sharper  interface  (/  =  0.1  nm)  is  also 
explored. 

5.3.  Simulations  and  results 

Considered  is  a  circular  twin  nucleus  of  initial  radius  ao  =  3  nm, 
embedded  in  a  much  larger  rectangular  domain  ^2  of  dimensions 
L  X  H.  Initially,  a  homogeneous  displacement  gradient  field  Vu 
is  applied  everywhere  in  ^2.  Particular  boundary  conditions  are 
applied  on  9^2,  as  discussed  in  more  detail  below.  The  order 
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Fig.  2.  Finite  element  mesh  and  initial  conditions,  undeformed  shape. 

parameter  and  local  displacements  are  then  relaxed  according  to 
the  conjugate  gradient  algorithm  discussed  in  Section  4,  as  a  local 
stationary  point  (specifically,  a  local  minimum)  of  the  total  free 
energy  of  the  system  is  sought.  If  the  magnitude  of  the  applied 
displacement  gradient  is  smaller  than  a  critical  value,  the  inclusion 
will  shrink  and  may  even  disappear,  the  latter  scenario  suggesting 
that  a  purely  homogeneous  elastic  response  corresponds  to  a 
minimum  energy  state.  If  the  magnitude  of  applied  displacement 
is  larger  than  a  critical  value,  the  inclusion  will  tend  to  grow  until 
it  reaches  a  final  size  limited  by  interaction  with  the  external 
boundary  of  the  domain.  Growth  does  not  necessarily  imply 
that  a  configuration  with  an  inclusion  of  finite  size  possesses 
less  total  energy  than  would  the  case  of  homogeneous  elastic 
deformation;  i.e.,  the  configuration  associated  with  an  inclusion 
can  be  metastable  and  may  constitute  a  relative  minimum  of  the 
total  energy. 

Recall  from  Fig.  1  that  0  denotes  the  orientation  of  the 
habit  plane,  i.e.,  the  plane  of  twinning  shear.  Specifically,  lattice 
orientation  vectors  are,  in  vector  form, 

s  =  [cos9,  sin^]^,  m=[-sin^,  cos^]^.  (69) 

Recall  also  from  (57)  that  the  mechanical  driving  force  for  twinning 
is  r  =  (7  :  (s  (g)  m).  From  (69),  it  follows  that  for  a  pure  shear  stress 
state  (e.g.,  cru  nonzero),  r  =  cru  (cos^  0 — sin^  ^) ;  for  a  pure  tension 
stress  state  (e.g.,  ^22  nonzero),  r  =  ^22  sin  0  cos  0. 

Three  kinds  of  boundary  conditions  are  considered.  The  first 
is  simple  shear  with  Dirichlet  conditions  on  the  order  parameter 
along 

{u^=yX2,U2  =  0,r]  =  0}  yXuX2edQ.  (70) 

The  second  is  shear  with  Neumann  (free)  conditions  on  traction 
and  order  parameter  along  the  lateral  edges  of  9^?: 

u,  =  yX2,  U2  =  0  V[Xi,(X2  =  0,H)], 
t  =  an  =  0,  n-Vr]  =  0  V[(Xi  =  0,  L),  X2]. 

The  third  is  simple  tension  with  Neumann  (free)  conditions  on 
traction  and  order  parameter  along  the  lateral  edges  of  9^?: 

Ui=0,  U2=SX2  V[Xi,(X2=0,H)], 
t  =  an  =  0,  n-Vr]  =  0  V[(Xi  =  0,  L),  X2]. 

In  (70)  and  (7 1 ),  y  is  the  magnitude  of  applied  shear,  equal  to  twice 
the  applied  shear  strain.  In  (72),  e  is  the  applied  stretch,  equal  to 
the  applied  tensile  strain. 

A  representative  finite  element  mesh  is  shown  in  Fig.  2.  The 
mesh  is  highly  refined  towards  the  center  where  the  inclusion 


is  initially  located.  In  simulations  that  follow,  numerous  elements 
resolve  the  thickness  of  the  boundary  of  the  inclusion  in 
its  equilibrium  configuration,  enabling  accurate  resolution  of 
spatial  gradients  of  the  order  parameter.  The  mesh  shown 
in  Fig.  2  includes  257,296  linear  triangular  elements.  Some 
simulations  were  repeated  with  a  finer  mesh  of  503,062  linear 
triangular  elements;  unless  noted  otherwise,  results  of  interest  are 
insensitive  to  mesh  density. 

The  initial  radius  of  the  twin  embryo  (inclusion)  is  set  at  Gq  = 
3  nm  for  the  following  reasons: 

•  According  to  the  analytical  sharp  interface  solution  (Ap¬ 
pendix  E),  a  bifurcation  from  circular  to  elliptical  shape  occurs 
for  a  radius  of  3. 17  nm,  corresponding  to  the  equality  conditions 
in  (E.9).  Therefore,  for  ease  of  validation  of  the  numerical  tech¬ 
nique  with  the  analytical  solution,  circular  inclusions  with  radii 
smaller  than  3.17  nm  must  be  considered.  Comparison  of  larger 
nuclei  would  require  consideration  of  elliptical  shapes  (i.e.,  as¬ 
pect  ratios  cd  <  1)  and  would  unduly  complicate  comparisons 
with  the  analytical  solution. 

•  According  to  atomic  simulations  of  twin  structures  and 
twinning  dislocations  in  Mg,  the  minimum  equilibrium  size 
of  a  (1011){1012}  twin  embryo  is  17  atomic  layers  (3.32  nm) 
thick  [13,14].  Shapes  of  twin  embryos  modeled  using  both  DFT 
and  empirical  MD  potentials  [13,14]  were  not  always  circular 
cylinders,  but  areas  of  twin  embryos  were  reasonably  close  to 
those  of  a  circle  of  radius  3  nm.  Unlike  atomic  simulations, 
the  present  phase  field  simulations  do  not  resolve  individual 
dislocation  lines,  but  length  scales  involved  are  about  the  same. 

•  The  equilibrium  thickness  I  of  the  twin  boundary  is  1  nm. 
If  ao  ^  3  nm  were  chosen,  the  boundary  thickness  would 
encroach  on  the  center  of  the  twin.  If  ao  ^  3  nm  were 
chosen,  computational  cost  associated  with  resolving  the  twin 
boundary  would  become  prohibitive. 

The  analytical  solution  in  (E.9)  with  inclusion  radius  =  02  = 
3  nm  corresponds  to  a  critical  stress  of  too/ iJt  = 
0.038.  This  solution  has  been  derived  for  an  inclusion  embedded 
in  an  infinite  domain,  with  a  sharp  interface  between  matrix 
and  inclusion.  Thus,  if  effects  of  finite  domain  size  and  a  diffuse 
interface  are  not  strong,  in  the  present  phase  field  simulations 
applied  shear  y  ^  0.038  should  result  in  growth  of  the  twin, 
while  y  0.038  should  result  in  decay  [for  cases  in  which  simple 
shear  boundary  conditions  (70)  are  applied].  Eor  shear  loading,  the 
applied  stress  is  defined  as  the  following  volume  average,  using  the 
divergence  theorem  and  the  first  of  (52): 

rA  =  ^  [  fiX2  dS  =  ^  [  cru  dQ.  (73) 

Shown  in  Eig.  3  are  simulation  results  for  applied  shear  (a) 
y  =  0.01  and  (b)  y  =  0.10.  In  each  case,  0  =  0  and  boundary 
conditions  (70)  apply.  Notice  that  decay  [Eig.  3(a),  y  =  0.01  < 
0.038]  to  a  homogeneous  medium  (no  twin)  and  growth  [Eig.  3(b), 
y  =  0.10  0.038]  occur  in  agreement  with  the  analytical 

solution.  Eor  y  =  0.01,  the  twin  nucleus  shrinks  and  disappears; 
for  y  =  0.10,  the  inclusion  grows  until  it  is  repelled  by  rigid  outer 
boundaries  9^2  where  r]  =  0is  prescribed. 

Next  sought  are  critical  applied  strains  for  twin  nucleation  for 
the  cases  listed  in  Table  2,  which  encompass  various  boundary 
conditions,  material  properties,  and  lattice  orientations.  Critical 
strains  are  obtained  in  practice  by  seeding  the  domain  with 
an  inclusion  of  initial  radius  gq  =  3  nm,  and  then  applying 
the  deformation  [y  for  shear,  e  for  tension  as  in  (70)-(72)]  in 
small  increments  of  magnitude  0.001.  Eor  each  increment,  system 
degrees  of  freedom  relax  in  conjunction  with  energy  minimization, 
performed  numerically  via  the  conjugate  gradient  technique.  It 
was  found  that  in  the  first  few  hundred  conjugate  gradient 
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(a)y  =  0.01. 


{b)y  =0.01. 


{c)y  =0.10.  {d)y  =  0.10. 


Fig.  3.  Order  parameter  p  and  shear  stress  <712:  (a),  (b)  y  =  0.01  (twin  disappears)  and  (c),  (d)  y  =  0.10  (twin  grows). 


Table  2 

Phase  field  simulations  and  predicted  critical  strain  and  stress  for  twin  nucleation. 
For  Case  4,  /cii/2  =  2k2  =  k.Fot  Case  6,  yc  =  y /2  and  ta  is  divided  by  2.  For  Case 
7,yc  =  s  and  ta  =  o  12,  with  o  the  average  normal  stress. 


Case 

e 

(rad) 

1 

(nm) 

/ril//C22 

l/H 

Boundary 

Loading 

Yc 

-Chill 

1 

0 

1 

1 

1 

Eq.  (70) 

shear 

0.050 

0.020 

2 

0 

1 

1 

1 

Eq.  (71) 

shear 

0.050 

0.000 

3 

0 

1 

1 

1.5 

Eq.  (70) 

shear 

0.049 

0.012 

4 

0 

1 

4 

1 

Eq.  (70) 

shear 

0.046 

0.023 

5 

0 

0.1 

1 

1 

Eq.  (70) 

shear 

0.049 

0.044 

6 

71  jO 

1 

1 

1 

Eq.  (70) 

shear 

0.050 

0.042 

7 

71  {A 

1 

1 

1 

Eq.  (72) 

tension 

0.045 

0.005 

iterations,  substantial  changes  in  size  and  shape  of  the  twin  occur 
in  response  to  the  prescription  of  an  initially  sharp  interface,  as  the 
interface  diffuses  to  a  width  commensurate  with  parameter  I  of  the 
phase  field  model.  In  subsequent  energy  minimization  iterations, 
the  twin  embryo  either  grows  or  decays,  with  magnitude  of  the 
rate  of  growth  or  decay  versus  number  of  iterations  increasing  with 
increasing  strain  difference  from  the  critical  strain.  The  minimum 
applied  strain  for  which  monotonic  decay  (to  null  size)  does 
not  occur  is  deemed  the  critical  strain  yc.  The  accuracy  of  yc 
obtained  via  this  method  is  limited  to  the  applied  strain  increment 
magnitude  of  0.001.  Qualitative  profiles  of  total  energy  versus 
inclusion  (twin)  size  for  two  values  of  applied  shear  are  shown  in 
Fig.  4.  At  a  given  shear  y,  when  a  <  ac  the  twin  will  tend  to  shrink 
and  disappear,  following  the  energy  surface  downhill  to  the  left. 
When  a  >  Gc  the  twin  will  tend  to  grow  until  impeded  by  external 
boundaries,  following  the  energy  surface  downhill  to  the  right.  For 
each  of  the  simulation  cases  considered  in  Table  2,  application  of 
shear  y  =  yc  —  0.001  resulted  in  shrinkage  and  disappearance  of 
the  initial  twin  nucleus  of  radius  3  nm. 


Fig.  4.  Sketch  of  system  energy  versus  inclusion  size  for  two  values  of  applied  shear. 

Critical  strains  yc  and  corresponding  average  shear  stresses  Ta 
are  listed  in  Table  2.  For  cases  1-5,  critical  stresses  are  defined 
as  volume  averages  Ta  in  (73)  at  the  applied  strain  corresponding 
to  growth.  For  case  6  (shear,  0  =  tt/O),  results  in  Table  2  are 
normalized  via  multiplication  by  the  factor  (cos^  0  —  sin^  =  1/2 
to  enable  comparison  with  cases  1-5.  For  case  7  (tension,  0  = 
7t  /4),  critical  strain  and  stress  are  defined  respectively  as  yc  =  ^ 
and  Xa  =  (sin  ^  cos  ^22  d^2  =  o  12. 

Fig.  5  shows  equilibrium  order  parameter  contours  for  each 
simulation  at  the  applied  critical  strain.  Fig.  6  shows  corresponding 
shear  stresses.  For  cases  1  -6,  the  local  shear  stress  shown  is  ori2 ;  for 
case  7,  the  stress  shown  is  that  acting  on  a  45°  plane,  (^22  —  ) /2. 

The  following  observations  from  Table  2  and  Figs.  5  and  6  are 
noteworthy: 

Case  1.  The  equilibrium  shape  of  the  twin  embryo  is  non-circular 
and  slightly  concave  as  a  result  of  the  applied  deformation  and 
boundary  conditions.  The  twin  grows  horizontally  until  repelled 


850 


J.D.  Clayton,  J.  Knap  /  Physica  D  240  (2011)  841-858 


Fig.  5.  Order  parameter  p  at  applied  critical  strains  corresponding  to  Table  2:  (a)  Case  1,  (b)  Case  2,  (c)  Case  3,  (d)  Case  4,  (e)  Case  5,  (f)  Case  6  and  (g)  Case  7. 
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Stress  [GPaJ 


0.89 

0.56 

0.23 


-0.10 


-0.43 


(c)  Case  3. 


(d)  Case  4. 


Stress  [GPaJ 


2.69 

1.76 

0.83 

-0.11 

-1.04 


(g)  Case  7. 


Fig.  6.  Shear  stress  at  applied  critical  strains  corresponding  to  Table  2:  (a)  Case  1,  (b)  Case  2,  (c)  Case  3,  (d)  Case  4,  (e)  Case  5,  (f)  Case  6  and  (g)  Case  7. 
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Fig.  7.  Profiles  of  order  parameter  along  horizontal  mid-plane  X2  =  H/2  and 
vertical  mid-plane  Xi  =  L/2  for  Case  1  (isotropic  surface  energy)  and  Case  4 
(anisotropic  surface  energy). 

by  the  lateral  boundaries  at  Xi  =  0,  L  The  critical  shear 

displacement  gradient  yc  =  0.05.  The  applied  average  shear  stress 
is  significantly  smaller  than  the  stress  that  would  be  obtained  for  a 

homogeneous  elastic  slab  with  no  inclusion,  i.e..rA/(/x)/c)  ~  0.4  < 

1.  Fig.  7  shows  the  value  of  the  order  parameter  along  reference 
coordinates  measured  from  the  center  of  the  inclusion:  X2  =  H/2 
(the  horizontal  mid-plane)  andXi  =  L/2  (the  vertical  mid-plane). 
From  Fig.  7,  along  either  mid-plane,  the  thickness  of  the  region  over 
which  OA  <  rj  <  0.9  is  very  close  to  parameter  I  =  1  nm. 

Case  2.  Here  the  lateral  boundaries  are  assigned  free  Neumann 
conditions.  The  twin  grows  horizontally  until  it  fully  spans  the 
domain.  The  applied  critical  strain  yc  =  0.05  is  identical  to  that 
of  Case  1,  but  the  average  stress  is  negligible:  tAlil^yc)  ~  0.  The 
transformation  strain  and  rotation  of  twin  fully  accommodate  the 
imposed  shear  deformation,  and  local  shear  stresses  vanish  except 
at  a  few  locations  near  twin  boundaries,  as  is  clear  from  Fig.  6(b). 

Case  3.  The  equilibrium  shape  of  the  twin  embryo  is  longer  in  the 
horizontal  (Xi )  direction,  but  otherwise  is  similar  to  that  of  Case  1. 
The  applied  critical  strain  yc  =  0.049  is  slightly  smaller  than  that 
of  Case  1  (see  Table  2),  and  the  ratio  of  average  stress  to  critical 
strain  is  ry\/(/xyc)  ~  0.24.  Increasing  the  aspect  ratio  L/H  of  the 
domain  from  1  to  1.5  apparently  has  a  small  effect  on  the  critical 
strain  for  homogeneous  nucleation. 

Case  4.  The  equilibrium  shape  of  the  twin  embryo  is  wider  (in  the 
Xi -direction)  and  flatter  (in  the  X2-direction)  than  that  of  Case  1,  as 
is  clear  from  profiles  of  the  order  parameter  in  Fig.  7.  The  thickness 
of  the  twin  boundary  is  larger  along  the  left  and  right  (shorter) 
sides  and  smaller  along  the  top  and  bottom  (longer)  sides,  as  also 
quantified  in  Fig.  7.  In  contrast  to  the  isotropic  cases  (e.g..  Case  1), 
prescription  of  anisotropic  surface  energy  (/cii/2  =  2k22  =  k) 
results  in  a  convex  rather  than  concave  shape.  The  applied  critical 
strain  yc  =  0.046  is  smaller  than  that  of  Case  1,  and  the  ratio  of 
average  stress  to  critical  strain  is  larger:  r^/C/xyc)  ~  0.5. 

Case  5.  The  equilibrium  shape  of  the  twin  embryo  is  that  of  a 
parallelogram,  confined  to  the  center  of  the  domain  where  the 
mesh  is  sufficiently  refined  to  resolve  the  twin  boundary  interface. 
The  thickness  of  the  interface  is  much  smaller  than  that  of  Case  1,  in 
agreement  with  parameter  choice  /  =  0.1  nm.  The  applied  critical 
strain  yc  =  0.049  is  very  close  to  that  of  Case  1,  and  the  average 
shear  stress  r^/C/xyc)  ~  0.9.  Shear  stress  concentrations  of  1. 5  GPa 
arise  at  four  corner  locations  along  the  twin-matrix  interface.  For 
this  case  only,  the  equilibrium  shape  is  mesh  sensitive.  The  critical 
strain  is  not  deemed  mesh  sensitive,  however,  since  applied  strains 
yc  <  0.049  resulted  in  decay  and  disappearance  of  the  initial  twin 
embryo. 

Case  6.  The  equilibrium  shape  of  the  twin  embryo  differs  from 
that  in  Case  1 ;  here,  the  twin  appears  rotated  such  that  one  axis. 


Fig.  8.  Interpolation  function  of  ( 9 )  [  35  ] . 

in  reference  coordinates,  is  aligned  normal  to  the  direction  s  of 
twinning  shear  (30°  from  horizontal).  The  applied  critical  strain 
yc  =  y  12  =  0.050  is  identical  to  that  of  Case  1,  and  Xhliiiyc)  ~ 
0.84.  The  applied  shear  stress  is  larger  in  Case  6  than  in  Case  1 
because  the  applied  shear  y  here  is  0.10  rather  than  0.05. 

Case  7.  The  equilibrium  shape  of  the  twin  embryo  differs  drastically 
from  that  in  Case  1.  The  twinned  region  is  large  and  nearly 
symmetric  about  vertical  mid-plane  Xi  =  L/2.  The  applied  critical 
tensile  strain  yc  =  s  =  0.045,  the  average  applied  stress  Ta  = 
0.005,  and  the  ratio  of  average  stress  to  critical  strain  tAlil^yc)  ~ 
0.11. 

In  results  discussed  thus  far,  interpolation  function  tp  of  (9)  [35] 
is  prescribed  as  or  =  3.  Exploratory  simulations  were  also 

conducted  with  O'  =  1  and  a  =  5;  corresponding  interpolation 
functions  are  plotted  in  Fig.  8.  Results  of  interest  are  not 
substantially  different  with  different  choices  of  a.  For  example,  as 
shown  in  Fig.  9(a),  profiles  of  the  order  parameter  versus  distance 
from  the  inclusion  along  mid-planes  Xi  =  L/2  and  X2  =  H/2  do 
not  depend  strongly  on  profiles  shift  only  slightly  to  the  left 
(i.e.,  towards  the  inclusion  center)  with  increasing  a.  The  middle 
value  explored,  a  =  3,  for  which  (9)  degenerates  to  a  cubic 
polynomial,  provides  the  anti-symmetry  conditions  (p{\  —  r])  = 
and  would  seem  most  physically  appropriate  for  modeling 
flat  interfaces.  For  a  =  1,3,  or  5,  profiles  of  transformation  shear 
(pyo  are  indistinguishable  in  Fig.  9(b). 

5.4.  Discussion 

Results  in  Section  5.3,  Table  2,  and  Figs.  5  and  6  suggest 
the  following  trends  regarding  nucleation  of  a  twin  embryo  of 
minimum  initial  radius  3  nm: 

•  The  critical  far-field  (i.e.,  applied)  shear  displacement  gradient 
for  nucleation  and  growth  is  0.045  <  yc  <  0.050,  and  appears 
fairly  insensitive  to  type  of  boundary  condition  (i.e.,  Dirichlet 
or  Neumann  in  (70)  or  (71)),  aspect  ratio  of  the  domain  ^2, 
surface  energy  anisotropy  /C11//C22,  twin  boundary  thickness  /, 
and  initial  lattice  orientation  0. 

•  For  shear  loading  conditions  perfectly  aligned  with  the 
geometry  of  the  twin  system  (0  =  0),  the  average  shear 
stress  after  nucleation  and  growth  is  0  <  Ta/jh  <  0.044,  is 
smaller  than  that  corresponding  to  what  would  be  observed 
for  homogeneous  elastic  shear  (i.e.,  for  no  twin),  and  can  be 
sensitive  to  type  of  boundary  condition. 

•  Equilibrium  and  near-equilibrium  shapes  of  twins  are  non¬ 
circular  and  are  sensitive  to  surface  energy  anisotropy  ratio 
X11//C22.  twin  boundary  thickness  parameter  /,  boundary 
conditions,  and  loading  mode  (e.g.,  shear  or  tension). 

Table  3  compares  twin  nucleation  criteria  resulting  from  the 
present  simulations  with  those  of  other  studies.  The  lower  bound 
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Fig.  9.  Profiles  of  (a)  order  parameter  p  and  (b)  twinning  shear  (pyo  along  vertical 
(Xi  =  L/2)  and  horizontal  (X2  =  H/2)  mid-planes  for  Case  1  (y  =  0.05)  for 
different  values  of  interpolation  parameter  a. 


Table  3 

Critical  shear  deformation  for  twin  nucleation  or  growth:  comparison  with  other 
models. 


Model 

Reference 

Remarks 

Yc 

Phase  field 

Elastic  inclusion 

This  work 
[30] 

Simulations  1-7  [  Table  2] 
Analytical  solution, 
infinite  medium 

0.045-0.050 

0.038 

Theoretical  stress 

Atomistic 

calculation 

[25,48] 

[23,24] 

Yc  =  yol{27T) 

Glide  of 

1/17[1011](1012) 
partial  dislocation 

0.021 

0.004 

on  the  present  results  is  in  fair  agreement  with  the  analytical 
model  of  Lee  and  Yoo  [30]  outlined  in  Appendix  E,  where  here 
the  notation  yc  =  tooZ/x.  Differences  between  phase  field 
results  and  the  analytical  elasticity  solution  of  Appendix  E  are 
apparently  due  to  resolution  of  finite  boundary  conditions  and 
finite  interfacial  thickness  in  the  former.  The  lower  bound  on 
critical  shear  strain  obtained  through  the  phase  field  approach 
exceeds  that  from  “theoretical  strength”  yol (lit)  associated  with 
motion  of  twinning  partial  dislocations  [25,48],  and  far  exceeds 
the  stress  from  empirical  lattice  statics  calculations  [23,24].  The 
latter  lattice  statics  approaches  consider  motion  of  a  single  partial 
dislocation  in  a  long  planar  twin  boundary  interface  and  do  not 
account  for  the  finite  size  and  curved  shape  of  the  twin  nucleus. 
This  may  account  for  the  discrepancies  in  values  between  atomic 
simulation  results  and  the  first  two  models  listed  in  Table  3. 

Eurther  qualitative  agreement  between  the  present  model  and 
analytical  elasticity  solutions  is  apparent  upon  consideration  of 
simulation  results  in  Eig.  3  and  the  sketch  in  Eig.  4.  Recall  that 
yc  corresponds  to  an  unstable  equilibrium  point,  i.e.,  a  saddle 
point  on  the  energy  surface  parameterized  in  terms  of  both 
size  and  shape  of  the  inclusion.  At  the  saddle  point,  the  total 
energy  is  minimized  with  respect  to  shape,  but  maximized  with 


respect  to  size  [31,32].  In  agreement  with  the  analytical  solution 
of  Appendix  E,  applied  strains  less  than  (or  greater  than)  the 
critical  strain  result  in  unstable  decay  (or  growth)  of  the  twin 
embryo.  The  saddle  point  configuration  is  not  attained  in  the 
conjugate  gradient  calculations  (witness  complete  disappearance 
or  substantial  growth  of  inclusions  in  Eig.  3)  because  the  algorithm 
seeks  energy  minima  rather  than  saddle  points.  Recall  that  the 
total  free  energy  is  a  non-convex  functional  of  the  order  parameter, 
as  is  obvious  from  the  contribution  of  the  double-well  potential 
function  in  (20).  It  is  expected  that  (28)  may  exhibit  one,  many,  or 
no  stationary  points  (i.e.,  minima,  maxima,  and/or  saddle  points), 
as  is  characteristic  of  phase  change  problems  involving  non-convex 
free  energy  functionals  [20,31].  Examination  of  the  rate  of  change 
of  total  free  energy  versus  iterations  for  each  simulation  indicate 
that,  after  a  sufficient  number  of  conjugate  gradient  iterations,  the 
state  thus  attained  corresponds  to  a  local,  and  possibly  global  free, 
energy  minimum.  The  same  end  result  could  presumably,  but  not 
necessarily,  be  obtained  via  use  of  a  kinetic  approach  to  evolve  the 
order  parameter  such  as  (D.3)  of  Appendix  D,  wherein  the  right 
side  of  (D.3)  vanishes  at  every  material  point  in  the  domain  in 
the  equilibrium  state  when  a  stationary  condition  on  the  total  free 
energy  has  been  attained. 

Remarks  regarding  assumptions  and  limitations  of  the  present 
phase  field  approach  are  in  order.  The  general  theory  developed 
in  Section  2  accounts  for  geometric  nonlinearity  (i.e.,  large 
deformations)  and  anisotropy.  Simulations  presented  in  Section  5 
incorporate  geometric  linearity  (Section  3)  and  isotropic  elastic 
constants.  Comparison  of  (39)  with  (57)  shows  that,  when 
elastic  strains  are  small,  the  difference  in  driving  force  for 
twinning  r  between  geometrically  nonlinear  and  linear  theories 
is  a  factor  of  which  in  turn  is  on  the  order  of  twinning 

shear  yo,  0.1295  for  the  twin  system  under  consideration  in  Mg 
(Table  1).  The  difference  between  the  two  theories  would  be 
more  severe  in  cubic  crystals,  where  for  the  usual  twinning  mode, 
yo  =  2“^/^.  As  remarked  in  Section  5.1,  anisotropy  of  elastic 
constants  in  Mg  is  low  [47].  Eurthermore,  atomic  simulations  have 
demonstrated  qualitative  and  reasonably  quantitative  agreement 
for  many  features  of  twins  in  hexagonal  crystals  modeled  using 
pair  potentials  [23]  which  cannot  correctly  describe  the  five 
independent  anisotropic  elastic  constants,  and  using  many-body 
potentials  [24]  which  can  correctly  describe  anisotropic  elastic 
constants.  The  present  phase  field  method  does  not  enable 
resolution  of  atomic  details  of  defect  structures  afforded  by 
quantum  or  molecular  mechanics  models.  However,  a  major 
advantage  of  the  present  approach  is  that  very  few  material 
parameters  are  needed:  at  minimum,  two  elastic  constants  and 
two  other  constants  related  to  twin  boundary  energy  and  twin 
boundary  thickness.  Eurthermore,  the  same  continuum  theory  can 
be  used  to  address  systems  of  much  larger  size  than  can  be  treated 
using  discrete  atomic  techniques,  so  long  as  the  numerical  grid  is 
sufficient  to  resolve  gradients  of  the  order  parameter  at  interfaces. 
The  present  work  has  considered  only  a  2D  idealization  of  twinning 
in  Mg.  Real  deformation  processes  can  involve  three  spatial 
dimensions  and  multiple  twin  systems,  up  to  six  (1011)  {1012} 
systems  in  Mg  [34].  Euture  work  may  consider  3D  simulations  and 
multiple  twins,  following  the  theory  in  Appendix  B. 

6.  Conclusions 

A  nonlinear  theory  has  been  developed  to  address  mechani¬ 
cal  twinning.  The  general  theory  accounts  for  large  deformations 
with  a  focus  on  equilibrium  thermodynamics.  A  geometrically  lin¬ 
ear  version  of  the  theory  has  been  implemented  in  finite  element 
simulations  of  homogeneous  twin  nucleation.  The  application  ap¬ 
pears  to  be  the  first  documented  phase  field  study  of  homogeneous 
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twin  nucleation  in  Mg  single  crystals.  Results  are  in  fair  agree¬ 
ment  with  analytical  solutions  for  critical  applied  strains  associ¬ 
ated  with  nucleation  and  subsequent  unstable  growth  of  a  small 
twin  embryo.  Critical  shear  displacement  gradients  for  nucleation 
and  growth  of  a  cylindrical  twin  of  size  3  nm  are  predicted  on 
the  order  of  0.045-0.050.  Numerical  simulations  enable  consider¬ 
ation  of  effects  of  boundary  conditions,  finite  domain  sizes,  and 
surface  energy  anisotropy  not  amenable  to  analytical  solutions. 
Results  demonstrate  that  minimum  applied  stresses  necessary  for 
twin  nucleation  and  growth  can  be  sensitive  to  boundary  condi¬ 
tions  associated  with  mechanical  and  order  parameter  fields,  and 
that  equilibrium  shapes  of  twin  nuclei  can  be  sensitive  to  surface 
energy  anisotropy. 


Letg  :  X  X  R"^^^  ^  R  be  defined 

g[X,  u,  Vu]  =  g[X,  (77,  X),  (V77,  F)].  (A.6) 

Assume  that  g  in  (A.6)  is  polyconvex  in  Vu.  Application  of 
definitions  (A.2)-(A.4),  with  n  =  3,  m  =  4,  nAm  =  3 
and  r(n,  m)  =  34,  yields  the  representation  g(X,  u,  Vu)  = 
g(X,  u,  Vu,  adj2  Vu,  adj3  Vu).  Array  adj2  Vu  contains  components 
of  cof  Vu^  and  determinants  of  all  2  x  2  sub-matrices  of  Vu  of  type 


Vi77  V277 
Fii  Fu 


(A.7) 


Similarly,  array  adj^Vu  contains  detF  and  determinants  of  the 
three  3x3  sub-matrices  of  Vu  of  type 


Appendix  A.  Existence  of  minimizers  of  energy  functional  (10) 

Existence  of  minimizers  of  energy  functionals  (10)  can  be 
demonstrated  under  conditions  of  convexity,  quasiconvexity  or 
polyconvexity  of  the  underlying  energy  density  functions.  The 
present  treatment  is  limited  to  convex  or  polyconvex  energy 
density  functions,  important  classes  of  functions  for  applications. 
Provided  here  is  a  brief  summary  of  essential  results;  textbooks  on 
the  calculus  of  variations  (e.g.,  [39])  should  be  consulted  for  a  more 
comprehensive  treatment. 


A 1.  Convex  energy  density  functions 

Let  W(X,  F,  77)  and  /(X,  77,  V77)  be  convex  in  F  and  V77, 
respectively.  In  addition,  let  g(X,  77,  F,  V77)  =  W(X,  F,  77)  -h 

/(X,  77,  V77).  It  is  easily  verified  that  g  is  a  convex  function  of 
(F,  V77).  Consequently,  it  can  be  shown  [38,39]  that  problem  (10) 
possesses  a  solution  provided  function  g  satisfies  the  following 
coercivity  condition: 

g(X,  t],  F,  Vt])  >  a(X)  +  Pi\¥\^  +  |Vj)|2)P/2  (A.1) 

for  almost  every  X  €  Q,  for  every  F  e  p  >  0  and  Vt]  e  M? 
and  for  some  absolutely  integrable  function  a  and  p  >  1.  |  •  | 
denotes  the  two-norm. 


A2.  Polyconvex  energy  density  functions 


A  function g  :  R"^^"  ^  R  =  RU{-hoo}  is  said  to  be  polyconvex 
if  there  exists  a  convex  function  h  :  ^  R,  such  that 

g(A)  =  h[T(A)],  (A.2) 

where  T  :  R"^^"  ^  ^jefined  as 

T(A)  =  (A,  adj2A,  . . . ,  adj„/,^A).  (A.3) 

Here,  n  A  m  =  min{n,  m},  adj^A  denotes  all  s  x  s  minors  of 
A  e  2  <  s  <  n  A  m,  and 


nAm 

r(n,  m)  =  ^ 

S=1 


mini 

(s!)2(m  -  s)!(n  -  s)l  ’ 


(A.4) 


In  the  context  of  conventional  strain  energy  density  functions 
W  :  R+^^  ^  R,  the  assumption  of  polyconvexity  yields  the  well- 
known  functional  form  W(F)  =  W(F,  cof  F^,  detF),  where  cof  A 
denotes  the  cofactor  of  matrix  A  [39,40],  and  where  W  is  a  convex 
function  of  its  arguments. 

Define  u=  (77,  x)  and  its  gradient 


Vi?y 

Vzt] 

'Vrj' 

fii 

Fu 

fl3 

F 

f21 

F22 

F23 

F31 

F32 

F33  _ 

"Vi77  V277  V377" 

Fll  Fi2  Fi3 

_  F21  F22  F23  _ 


(A.8) 


It  is  emphasized  thatg  is  convex  in  its  last  three  arguments  for  all 
Xe  C2  and  u  g  R"^  by  the  poly  convexity  ofg. 

Notions  of  polyconvexity  become  further  simplified  if  an 
additive  decomposition  of  g  is  used.  Take  W  polyconvex  in  F, 
W(X,  F,  77)  =  W(X,  F,  cof  F^,  detF,  77),  and  /  convex  in  V77; 

function  g  is  then  represented  as 

g(X,  u,  Vu)  =  g(X,  77,  F,  cof  F^,  detF,  V77) 

=  W(X,  F,  cof F^,  detF,  77)  -h/(X,  77,  V77),  (A.9) 

and  clearly  satisfies  the  more  general  polyconvex  form  above. 
Furthermore,  since  W  and  /  are  both  convex  functions  of  their 
appropriate  arguments,  g  is  a  convex  function  of  F,  cof  F^,  detF, 
and  V77,  which,  in  turn,  renders  g  polyconvex. 

The  existence  of  minimizers  of  (10)  with  g  in  (A.9)  is  assured 
provided  the  following  coercivity  condition  is  satisfied  [39]: 


3 

gix,  u,  Vu)  >  a(X)  +  Aladj,Vu|P^  (A.10) 

S=1 

where  a  is  an  absolutely  integrable  function  over  C2,  Ps  >  0. 
Pi  >  2,  P2  >  and  p3  >  1.  This  condition  can  be  expressed 
more  explicitly  as 

giX,  u,  Vu)  >  a(X)  +  P,(\F\^  +  |V/7p)P'/2 

+  ft|cofF|P2 +^3(detF)P=’.  (A.11) 


Appendix  B.  Multiple  twin  systems 

The  theory  of  Section  2  is  extended  to  account  for  multiple  twin 
systems,  i.e.,  more  than  two  phases.  Denote  scalar  order  parameter 
fields  by  rjifX,  t).  Let  i  =  0,  1,  2, . . . ,  n,  with  n  -h  1  the  total 
number  of  phases.  For  a  crystal  with  multiple  twin  systems,  n  is 
the  number  of  twin  systems,  i  =  0  is  associated  with  the  parent, 
and  i  =  1,  2, . . . ,  n  are  associated  with  twin  systems.  Denote 

770  =  1  VX  G  parent,  77/  =  1  VX  g  twin  system  i, 

0  <  77/  <  1  VX  G  interfaces.  (B.l) 

The  following  constraint  applies  [3]: 

n 

^?jj(X,f)  =  l.  (B.2) 

i=0 

Order  parameter  77/  can  be  interpreted  as  the  mass  fraction  of  phase 
i  at  a  given  material  point  X.  Twinning  kinematics  for  a  given 
system  i  are  described  by  the  simple  shear 

Fj  =  1  -h  yojSi  (g)  m/  (with  S/  •  m/  =  0,  S/  •  S/  =  m/  •  m/  =  1). 

(B.3) 


[Vu]  = 


(A.5) 
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The  unit  normal  to  the  surface  of  composition  (i.e.,  the  habit  plane) 
in  the  reference  configuration  is  m/.  The  magnitude  of  the  twinning 
shear  and  the  shear  direction  are  /o/  ^nd  Sj,  respectively.  The 
deformation  gradient  F  =  V/  is  decomposed  as  in  (6): 

F  =  F^F^  (B.4) 

where 

F^(X,  t)  =  F(F'^)“^  =  elastic  deformation, 

F^[r]i(X,  t);  i  =  1,  2, . . . ,  n]  =  stress-free  twinning  shear. 
Specifically,  twinning  shear  is  interpolated  as  follows: 

n 

F''  =  1  +  '^[<Pi(}?i)]roiSi  <8  nii,  (B.6) 

2=1 

where  interpolation  functions  analogous  to  (9)  are 

(Pi(rii)  =  ar]f  -h  2(2  -  a)r]f  -h  (o'  -  3)rif,  0  <  a  <  6.  (B.7) 

Clearly  ¥^(rii  =  1;  i  >  0)  =  1  -h  yo/Si  0  m/  =  F,  and  =  OVi  > 
0)  =  1.  For  a  single  twin  system,  det  F^  =  det  (1  -h  cpiYoiSi  <S>  m/)  = 
1  +  (PiYoi^i  •  m/  =  1,  so  volume  is  conserved.  However, 

—  (detF'')  =  (detF’>)tr[F''(F'')“’] 
df 

n 

^  (det  ¥’>)  ^  (pYoi  Si  •  m;  =  0,  (B.8) 

2=1 

SO  volume  is  only  conserved  approximately  by  twinning  when 
multiple  phases  are  present  (e.g.,  at  triple  points).  Strain  energy 
per  unit  reference  volume  W  is  of  the  functional  form 

W  =  t^i),  (B.9) 

where  is  the  defined  in  (11).  Strain  energy  density  W  always 
vanishes  at  null  elastic  strain: 


W(0,  •)  =  0.  (B.IO) 

Strain  energy  density  and  second-order  moduli  are  written 


E^=0 


W  =  -E^  :  C(.0  :  E^  C(,0  = 

For  a  compound  twin  in  a  centrosymmetric  lattice, 

Qj  =  2mi  (g)  nif  —  1, 


(B.ll) 

(B.12) 


and  elastic  coefficients  of  twins  C/  are  related  to  those  of  the  parent 
Co  by 


Qabcd  —  —  V)abcd 

=  OiAE  QiBF  OiCG  OiDH  QeFGH  • 


(B.13) 


Elastic  coefficients  in  the  interface  are  interpolated  the  same  way 
as  the  twinning  shear: 


n 

C(r]i)  =  Co  +  ^[Ci  -  Co]q>i.  (B.14) 

2=1 

In  the  isotropic  approximation,  Co  =  Cj  Vi  and  W  does  not  depend 
explicitly  on  The  local  interfacial  energy  per  unit  reference 
volume  is  written  as  follows,  extending  (18): 

n 

fivu  V??j)  =  +  iCi  ■  ®  V?;,)] 

2=1 


n-1  n 

gijitli,  rij),  (B.15) 

i=l  j=i+l 

with  Ki  symmetric  second-order  tensors  of  material  constants.  A 
possibility  for/oj  is 

/o,•=/^^)f(l-/7,■)^  (B.16) 


with  A  a  non-negative  constant.  The  following  function  penalizes 
more  than  two  phases  at  a  boundary  (e.g.,  at  triple  points)  [20]: 

gii(.r)u  r]j)  =  -  r),  -  rjj),  (B.17) 

with  B  a  non-negative  constant.  Energy  function  (B.15)  degenerates 
to  (18)  when  only  one  twin  system  is  present.  Other  possible 
functions  that  might  be  considered  are  given  by  Steinbach  et  al.  [3]. 
The  total  free  energy  functional  ^  is  defined  as  the  volume  integral 

<!'=  WiE^,r]i)<iQ+  /  /(??,-,  V??,)  (B.18) 

J  J  C2 

Applying  Hamilton’s  principle,  the  following  variational  equation 
is  posited  to  suggest  static  equilibrium  and  boundary  conditions: 


8^  - 


t  •  d5 


hi8r}i  d5  =  0, 


(B.19) 


where  hi  is  a  scalar  conjugate  force  to  variation  in  order  parameter 
Tji.  Taking  the  first  variation  of  the  interfacial  energy  and  applying 
the  divergence  theorem,  the  analog  of  (29)  is 


Ja  i^  dVi 


-  (  2k:,  :  [V (V r]i)]Sr]i  di2  +  f  2k:,  :  [(V?j,)  (g)  n]  Sr],  dS 
Js2  JdQ  \ 


n— 1  n  p 

+  /  -^Sr]id^2. 


Bgij , 


(B.20) 


Taking  the  variation  of  strain  energy  with  F  and  rji  the  independent 
variables. 


8  Wd^ 
Jq 


aw 

—8r]i  d^2 
drii 


dW~ 


•  ax  d^2 


axd5. 


(B.21) 


For  admissible  variations  ^x  and  8r]i,  it  follows  that  in  ^2, 
Euler-Lagrange  equations  are 


V 


aw 


=  V  ■P  =  0, 

m 


dr]i 


2iCi  :  [V(V/72)]  ■ 


aw 

drji 


-E 

j=i+l 


dgij 

drii 


(i  =  1,  2, . . . ,  n). 

Boundary  conditions  are 

t  =  Pn;  hi  =  2iCi  :  (Vr]i  0  n)  (i  =  1,  2,  . . . ,  n). 


(B.22) 

(B.23) 

(B.24) 


Appendix  C.  Interfacial  energy 


Consider  the  theory  of  Sections  2  and  3  wherein  a  single  order 
parameter  is  sufficient.  Twinned  region  and  parent  are  separated 
by  a  thin  interface  in  which  0  <  rj  <  1.  The  contribution  of  the 
phase  field  to  the  free  energy  density,  in  the  absence  of  elastic 
strain  energy,  is  written  as  the  Taylor  series  expansion  [  1  ] 


f(r],  Vt],  VVt],  ...)  =  foil))  -h 


1 

+  - 


9/ 

dVrj 

d^f 


9/ 

•  Vyy  H - 

0  '  avv^ 


:  VV?7 


2  dVrjdVrj 


(Vrj  (g)  Vrj)  -h  . 


■■fo(v)  +  A-  Vr]  +  ^:  (VVt]) 


+  ^2;:  (V/j®  V/7)  +  ...,  (Cl) 
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where  zero  subscripts  correspond  to  a  uniform  value  of  the  order 
parameter.  Here  it  is  assumed  for  convenience  that  minima  at 
uniform  composition  occur  at  wells /o(0)  =  /o(l)  =  0.  and  that 
/o  >  0;  in  terminology  of  Allen  and  Cahn  [2],  A/o  =  /o(«)  — /o(0)  = 
/o(»)  — /o(l)  =  /o(*)-  The  analysis  is  limited  to  centrosymmetric 
structures  for  which 


A  =  0. 


(C.2) 


Other  coefficients  can  depend  on  the  order  parameter  [1].  For 
isotropic  or  cubic  crystals,  §  =  and  ^  f  1.  Integrating  (Cl) 
over  ^2  and  omitting  higher-order  terms, 


[  fdQ 

JC2 


!.[ 


/o  +  I :  (VVj))  +  :  (V/j  (gi  V?)) 


AS2. 


(C.3) 


Applying  the  divergence  theorem  and  the  chain  rule, 
I  §  :  (VV^)  dQ  =  I  §  :  (Vt]  (8)  n)  d5 

Jc2  JdQ 

-  f  (V  •  §)  •  (V/y)  d^2 


=-/ 


drj 


:  (V?7  (g)  V?7)  di7, 


(C.4) 


assuming  Vr]  =  0  along  9^2  [1].  Substituting  (C4)  into  (C3), 


\fo  + ic  :  (Vr]  ^Vr])]d^2, 
J  Q 

where 


(C.5) 


K  = 


drj 


i  1  9Y 

92/ 

'  0  2  dVrjdVrj 

0  drjdVVr]  0 

(C6) 


For  cubic  or  isotropic  symmetry,  (C.5)  and  (C.6)  reduce  to  [1] 


[  [fo  +  K\Vr]\^]dn, 

Jr2 


K  = 


fi. 

FI] 

1  d^f 

92/ 

\2 

dr]J 

0  2  9|Vr,|2 

Q  drjdV^T] 

(C.7) 


Consider  a  flat  interface,  for  which  are  sought  the  interfacial  energy 
and  equilibrium  thickness  under  stress-free  conditions  [1,2].  In  the 
present  context,  these  are  physically  related,  respectively,  to  twin 
boundary  energy  per  unit  area  and  the  thickness  of  the  region 
near  the  boundary  where  atomic  positions  deviate  from  those  of 
a  perfect  lattice.  Henceforth  in  Appendix  C,  it  is  assumed  that  k  = 
k\  =  constant.  At  equilibrium,  when  strain  energy  W  vanishes, 
(32)  gives 


-f  =  iKV^r^.  (C.8) 

dr] 

Consider  an  infinitely  long  columnar  volume  with  element 
d^2  =  AodX  having  constant  cross-sectional  area  Aq  and  normal 
Lagrangian  coordinate  to  the  interface  X.  Assume  that  r]  varies  only 
in  theX  direction:  rj  =  ri(X).  Then  (C.8)  becomes 


^=2k^ 

dr]  9X2  • 

Integrating  both  sides  and  applying  the  chain  rule  [2], 


Jo  dr] 


{d[Kidr]/dXy]/dr]}dr] 


^  dr]/dX  =  (fo/tcy^\ 


(C.9) 


(CIO) 


Interfacial  energy  per  unit  area  F  is  defined  as  [1] 

/+00  p+oo 

fdX=  \fo  +  K(dr]/dXy]dX 

-oo  J  — oo 

/+00  p  1 

/o  dx  =  2  /  (/c/o)’/2  dr], 

-OO  J  0 


(C.11) 


where  limits  ri(X  — oo)  =  0  and  ri(X  +oo)  =  1.  A 

central  difference  approximation  of  the  equilibrium  thickness  I  of 
the  interfacial  region  is  [  1  ] 

dn 

^  ^  [r](X+)  -  r](X-)]/[X+  -  X-]  «  1/1 


«  (/o,max//r)’/^  (C.12) 

where  for  a  symmetric  potential /o  with  maximum  at  the  midpoint 
of  the  interface. 


/o,max  —fo\x-+l/2  —fo\x+-l/2  — /ol?7=l/2-  (C.13) 

Applying  this  argument  to  the  double-well  function/o  of  (20)  gives 
the  “height”  of  the  well: 

/o, max  =/o(l/2)=  A/16.  (C.14) 

From  (C.12),  the  equilibrium  thickness  of  the  interface  is 


(  «  (^//o.max)’^^  =  (C.15) 

From  (C.ll),  the  equilibrium  surface  energy  of  the  interface  is 

(/c/o)’/2  d/?  =  2(a/c)’/2  l\(\-r])dr] 


r 


T 


iAKy^^/3. 


Use  of  (C.15)  and  (C.16)  results  in  k 
listed  in  (21). 


3r//4  and  A 


(C.16) 

12r//as 


Appendix  D.  Kinetic  approach 


A  typical  approach  for  modeling  time  evolution  of  the  order 
parameter  towards  an  equilibrium  point  involves  replacement  of 
static  equilibrium  condition  (32)  with  a  kinetic  equation,  following 
the  TDGL  formalism.  Phase  equilibrium  can  be  expressed  as 

8i//8r]  =  dfo/drj  -  2ic  :  [V(V^)]  -h  (dW/dr])\f  =  0.  (D.l) 

The  kinetic  approach  [2]  suggests  an  evolution  law  of  the  form 
(drj/dt)\x  =  rj  =  -l[8ir/8rj}\  L  =  constant  >  0.  (D.2) 

Substituting  gives 

ij  =  -L{dfo/drj  -  2k  :  [V(Vrj)]  +  (dW/drj)\,}.  (D.3) 

This  evolution  equation  is  supplemented  by  initial  conditions 
r](X,  0)  throughout  domain  ^2.  Equilibrium  is  approached  as  rj  ^ 
0  [2].  A  similar  approach  can  be  applied  in  the  context  of  multiple 
order  parameters  (Appendix  B,  [3,20]).  In  that  case,  from  (B.23), 
evolution  laws  of  the  following  form  are  suggested,  where  4/  =  Uk 
is  a  positive  definite  matrix  of  kinetic  coefficients: 


drji 


2iCi  :  [V(Vrjd] 


dr]i  F  dr]i 


(D.4) 


Appendix  E.  Analytical  solution,  homogeneous  twin  nucleation 

Following  Lee  and  Yoo  [30],  considered  in  what  follows  are 
conditions  for  nucleation  of  a  two-dimensional  (i.e.,  cylindrical) 
twin  embedded  in  an  otherwise  homogeneous  elastic  solid  of 
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infinite  extent.  The  analysis  of  Lee  and  Yoo  [30]  builds  on  earlier 
work  of  Eshelby  [42]  and  Johnson  and  Cahn  [49].  Similar  analyses 
of  the  Gibbs  free  energy  change  associated  with  twin  nucleation 
have  since  appeared  [11.31-33].  The  present  analysis  is  limited  to 
isotropic  linear  elasticity,  with  identical  elastic  constants  shared 
by  the  twin  (i.e.,  the  inclusion)  and  the  parent  (i.e.,  the  matrix). 
Consider  an  ellipsoidal  twin  with  semi-axes  ai,  02,  and  as,  where 
Os  ^  00  such  that  the  ellipsoid  degenerates  to  an  infinitely 
extended  cylinder.  The  cross-sectional  area  of  the  twin  is  At  = 
7ta^a2,  and  the  aspect  ratio  of  the  twin  is  co  =  02/0^  <  1.  The  total 
Gibbs  free  energy  change  per  unit  length  of  the  cylinder  associated 
with  twin  nucleation  is 


AG  =  We -h  ~ ‘^t'^ooVo-  (E-^) 

The  elastic  strain  energy  per  unit  length  is  [41] 


We  =ATS(o/(-l+(oy,  S  =  At/oV(2  -  2v).  (E.2) 

The  surface  energy  per  unit  length  of  the  twin  embryo  is  [30] 

0s  =  4r[AT/(7tcv)]E(l<,  7r/2),  (E.3) 

where  the  complete  elliptic  integral  of  the  second  kind  with  k  = 
(1  —  is 


(E.4) 


The  work  per  unit  length  associated  with  far-field  loading  stress 
(Too  is  [30.31] 


—  AtTooYo  =  —AtIcToo  •  (s0in)]yo-  (E.5) 

The  critical  aspect  ratio  and  critical  size  for  twin  nucleation  for 
a  given  far-field  stress  are  determined  by  the  simultaneous 
stationary  conditions 


dAG 


=  0, 


dAG 


=  0. 


do)  dAj 

Using  (E.1)-(E.5),  these  two  conditions  can  be  expressed  as  [30] 


(E.6) 


^oo/oA  =  co/(l coy co(co  -  !)£ 

X  {(1  -h  coy[2co(dE/dco)  -  E]}-\  (E.7) 


At  =  {2r(l-E(oy[2(o(dE/d(o) -E] 

Simultaneous  solution  of  (E.7)  and  (E.8)  gives  the  critical  aspect 
ratiocaandsize,  i.e.,  area^T  or  major  semi-axis  ai  =  [^^/(Tra;)]^/^, 
of  a  twin  nucleus  for  a  given  far-field  stress  too  -  Noting  the  solution 
is  a  saddle  point  [31,32]  described  by  d^(AG)/d(o^  >  0  and 
d^(AG)/dAj  <  0.  it  follows  that  the  solution  is  unstable  with 
respect  to  changes  in  size  of  the  inclusion.  Eor  a  fixed  aspect  ratio 
and  fixed  far-field  stress,  if  an  embryo  smaller  than  the  critical  size 
given  by  (E.8)  is  introduced  in  the  matrix,  it  will  tend  to  disappear, 
while  if  it  is  larger  than  the  critical  size,  it  will  tend  to  grow.  The 
greater  the  applied  stress  too,  the  smaller  the  critical  size,  meaning 
that  given  a  statistical  distribution  of  potential  twin  embryo  sizes 
and  shapes,  nucleation  becomes  more  feasible  as  the  applied  stress 
increases  [31,32].  When  too  exceeds  the  following  threshold,  the 
solution  is  a  circular  cylinder  [30]: 


too  >  55/(12yo)  =  5/x]/o/[24(l  -  y)]  :  a;  =  1, 
At  =  Ttr^lTooVo  -  S/4]~\ 
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